#### FIGURE C.1: POLICY RIGHTS INDEX
#### Pooled effects on policy rights index

rm(list = ls())
source("./2_code/00_setup.R")

#### STUDY 1 ####

data <- fread(paste0(data_path, "data_study1.csv"), header = TRUE)
data$worked <- as.factor(data$worked)

# Narratives on policy index
lin_policy_economic <- lm_lin(policies_simple_s ~ treat_economic, covariates = ~
                                edad + worked + male, data = data[(treat=='economic' | treat=='control')])
lin_policy_humanitarian <- lm_lin(policies_simple_s ~ treat_humanitarian, covariates = ~
                                    edad + worked + male, data = data[(treat=='humanitarian' | treat=='control')])


#### STUDY 2 ####

data2 <- fread(paste0(data_path, "data_study2.csv"), header = TRUE)
data2$worked <- as.factor(data2$worked)

# Narratives on policy index
lin_policy_economic2 <- lm_lin(policies_simple_s ~ treat_economic, covariates = ~
                                 edad + worked + male, data = data2[(treat=='economic' | treat=='control')])
lin_policy_humanitarian2 <- lm_lin(policies_simple_s ~ treat_humanitarian, covariates = ~
                                     edad + worked + male, data = data2[(treat=='humanitarian' | treat=='control')])


#### STUDY 3 ####

data3 <- fread(paste0(data_path, "data_study3.csv"), header = TRUE)
data3$employed <- as.factor(data3$employed)

# Narratives on policy index
lin_policy_economic3 <- lm_lin(policies_simple_s ~ treat_economic, covariates = ~
                                 edad + employed + male, data = data3[(treat=='economic' | treat=='control')])
lin_policy_humanitarian3 <- lm_lin(policies_simple_s ~ treat_humanitarian, covariates = ~
                                     edad + employed + male, data = data3[(treat=='humanitarian' | treat=='control')])


#### META ANALYSIS ####

### Calculate observations for economic treatment ###

dataset_names <- c("data", "data2", "data3")

calculate_total_observations_eco <- function(dataset) {
  dataset %>%
    filter(treat == 'economic' | treat == 'control') %>%
    summarise(total_observations = sum(complete.cases(across(c(policies_simple_s)))))
}

total_observations <- lapply(dataset_names, function(name) {
  calculate_total_observations_eco(get(name))
})

obs_eco <- unlist(total_observations)

types_economic <- data.frame(
  treatment = rep('economic', 3),
  Study = c('Study 1', 'Study 2', 'Study 3'),
  n_total = obs_eco
)


### Calculate observations for humanitarian treatment ###

calculate_total_observations_hum <- function(dataset) {
  dataset %>%
    filter(treat == 'humanitarian' | treat == 'control') %>%
    summarise(total_observations = sum(complete.cases(across(c(policies_simple_s)))))
}

total_observations_hum <- lapply(dataset_names, function(name) {
  calculate_total_observations_hum(get(name))
})

obs_hum <- unlist(total_observations_hum)

types_humanitarian <- data.frame(
  treatment = rep('humanitarian', 3),
  Study = c('Study 1', 'Study 2', 'Study 3'),
  n_total = obs_hum
)


### Economic treatment, policy index outcome ###

out_economic <- rbind(
  tidy(lin_policy_economic) %>% filter(term=='treat_economic') %>% dplyr::select(estimate, std.error),
  tidy(lin_policy_economic2) %>% filter(term=='treat_economic') %>% dplyr::select(estimate, std.error),
  tidy(lin_policy_economic3) %>% filter(term=='treat_economic') %>% dplyr::select(estimate, std.error)
)

out_economic <- cbind(out_economic, types_economic)

m.gen_eco <- metagen(
  TE = estimate,
  seTE = std.error,
  studlab = Study,
  data = out_economic,
  common = FALSE,
  random = TRUE,
  method.tau = "REML",
  n.e = n_total,
  title = "economic_treatment_policy"
)

summary(m.gen_eco)


### Humanitarian treatment, policy index outcome ###

out_humanitarian <- rbind(
  tidy(lin_policy_humanitarian) %>% filter(term=='treat_humanitarian') %>% dplyr::select(estimate, std.error),
  tidy(lin_policy_humanitarian2) %>% filter(term=='treat_humanitarian') %>% dplyr::select(estimate, std.error),
  tidy(lin_policy_humanitarian3) %>% filter(term=='treat_humanitarian') %>% dplyr::select(estimate, std.error)
)

out_humanitarian <- cbind(out_humanitarian, types_humanitarian)

m.gen_hum <- metagen(
  TE = estimate,
  seTE = std.error,
  studlab = Study,
  data = out_humanitarian,
  common = FALSE,
  random = TRUE,
  method.tau = "REML",
  n.e = n_total,
  title = "humanitarian_treatment_policy"
)

summary(m.gen_hum)


#### PLOTTING ####

# Extract data for economic treatment, policy index outcome
total_observations <- sum(m.gen_eco$n.e)
forest_data_eco <- data.frame(
  Study = c(m.gen_eco$studlab, "Pooled Effect"),
  Coefficient = c(m.gen_eco$TE, m.gen_eco$TE.random),
  Observations = c(m.gen_eco$n.e, total_observations),
  SE = c(m.gen_eco$seTE, m.gen_eco$seTE.random),
  Lower95 = c(m.gen_eco$lower, m.gen_eco$lower.random),
  Upper95 = c(m.gen_eco$upper, m.gen_eco$upper.random),
  Lower90 = c(m.gen_eco$TE - (1.645 * m.gen_eco$seTE), m.gen_eco$TE.random - (1.645 * m.gen_eco$seTE.random)),
  Upper90 = c(m.gen_eco$TE + (1.645 * m.gen_eco$seTE), m.gen_eco$TE.random + (1.645 * m.gen_eco$seTE.random)),
  Narrative = "Economic Narrative",
  Outcome = "Policy"
)

# Extract data for humanitarian treatment, policy index outcome
total_observations <- sum(m.gen_hum$n.e)
forest_data_hum <- data.frame(
  Study = c(m.gen_hum$studlab, "Pooled Effect"),
  Coefficient = c(m.gen_hum$TE, m.gen_hum$TE.random),
  Observations = c(m.gen_hum$n.e, total_observations),
  SE = c(m.gen_hum$seTE, m.gen_hum$seTE.random),
  Lower95 = c(m.gen_hum$lower, m.gen_hum$lower.random),
  Upper95 = c(m.gen_hum$upper, m.gen_hum$upper.random),
  Lower90 = c(m.gen_hum$TE - (1.645 * m.gen_hum$seTE), m.gen_hum$TE.random - (1.645 * m.gen_hum$seTE.random)),
  Upper90 = c(m.gen_hum$TE + (1.645 * m.gen_hum$seTE), m.gen_hum$TE.random + (1.645 * m.gen_hum$seTE.random)),
  Narrative = "Humanitarian Narrative",
  Outcome = "Policy"
)


### Combined plot ###

forest_data_combined <- rbind(forest_data_eco, forest_data_hum)

desired_order <- c("Pooled Effect", "Study 3", "Study 2", "Study 1")

forest_data_combined <- forest_data_combined %>%
  mutate(Study = factor(Study, levels = desired_order))

policy_plot <- ggplot(forest_data_combined, aes(x = Coefficient, y = Study)) +
  geom_point(aes(color = Study), size = 3) +
  geom_errorbarh(aes(xmin = Lower95, xmax = Upper95, color = Study), height = 0, linewidth = 0.9) +
  geom_errorbarh(aes(xmin = Lower90, xmax = Upper90, color = Study), height = 0, linewidth = 1.9) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "red") +
  labs(title = "", x = "", y = "", caption = "ATE Policy Index") +
  theme_minimal() +
  theme(
    axis.title.x = element_blank(),
    plot.caption = element_text(size = 11, hjust = 0.5),
    strip.text = element_text(size = 16),
    axis.text.y = element_text(size = 13),
    legend.position = "none"
  ) +
  scale_x_continuous(limits = c(-0.1, 0.3)) +
  scale_color_manual(values = c("Study 1" = "grey50", "Study 2" = "grey50", "Study 3" = "grey50", "Pooled Effect" = "red")) +
  facet_grid(~Narrative)

policy_plot

ggsave(filename = paste0(plot_path, "figure_C1.png"), plot = policy_plot, width = 8, height = 5, dpi = 600)
